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МЕТОД ГЕНЕТИЧЕСКОГО ПРОГРАММИРОВАНИЯ С СЕТЕВЫМ ОПЕРАТОРОМ 
ДЛЯ ИДЕНТИФИКАЦИИ СИСТЕМ УПРАВЛЕНИЯ 


Рассмотрена задача структурно-параметрической идентификации системы управления. Приведено ее решение 
по методу генетического программирования на основе сетевого оператора. Описаны свойства сетевого 
оператора, правила его построения и использования. Изложен принцип базисного решения, который позволяет 
применять генетические операции к вариациям такого решения. Показано влияние выбора базисного 
решения, задаваемого разработчиком на основе опыта исследования конкретной задачи, на эффек ТИВНОСТЬ 
генетического алгоритма. Приведены результаты вычислительного эксперимента решения прикладной 
задачи идентификации. 
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Введение. Решение задачи структурно-параметрической идентификации системы управления, 
как правило, проводится в два этапа: первоначально выбирается структура системы управления, 
а затем — значения параметров в ней. Этап выбора структуры в большинстве случаев 
осуществляется «вручную», т. е. в результате анализа причин наблюдаемого поведения системы. 
Вычислительные методы применяются только на этапе выбора значений параметров. 

Метод генетического программирования [1], разработанный в 1992 г. профессором 
университета Стэнфорда Дж. Козой, позволяет использовать вычислительную машину для поиска 
алгоритма или математического выражения, т. е. с помощью вычислительной машины определить 
структуру системы управления. Основным элементом генетического программирования является 
структура данных — польская запись, — которая позволяет описывать в виде строки символов 
любые математические выражения или алгоритмы. Генетический алгоритм применяет 
генетические операции, прежде всего операцию скрещивания, к польским записям и получает 
новые записи, которые соответствуют новым математическим выражениям. 

Для корректного использования генетического программирования требуется анализ 
достаточно большого множества возможных решений. При проверке каждого решения, т. е. при 
вычислении функции приспособленности для него или при выполнении генетических операций с 
ним, каждая польская запись должна быть корректно прочитана, для чего должен использоваться 
лексический анализатор. Работа лексического анализатора замедляет процесс вычислений. В 
Стэнфорде при решении задач методом генетического программирования используют мощную 
вычислительную машину из тысячи кластеров [5]. Другим недостатком генетического 
программирования, основанного на польской записи, является то, что сама запись первоначально 
не предусматривает включения в нее несколько раз одних и тех же аргументов, переменных или 
констант. Наконец, третьим недостатком генетического программирования с польской записью 
является то, что при ее использовании трудно ограничить область поиска. В примерах 
использования генетического программирования [1—4] размеры первоначальных популяций 
превосходили один миллион. 

В генетическом программировании альтернативой польской записи является сетевой 
оператор. Это другая структура данных, которая позволяет описывать математические выражения 
с помощью ориентированного графа [6—10]. Данный граф может представляться в 
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вычислительной машине в виде целочисленной верхней треугольной матрицы. Сетевой оператор 
лишен перечисленных выше недостатков польской записи. Поэтому, если исследователю 
предварительно известен приблизительный вид искомого математического выражения, то он 
может задать некоторое базисное решение и выполнить поиск в его окрестности, что значительно 
сужает область поиска и уменьшает размер первоначальной популяции. 
Постановка задачи. Рассмотрим систему автоматического управления, которая состоит из двух 
подсистем. Математическая модель одной из них известна. Математическую модель второй 
подсистемы требуется определить. Известны размерности каждой из подсистем и то, что модель 
неизвестной подсистемы может быть описана обыкновенными дифференциальными уравнениями. 
Необходимо по вектору управления, действующему на систему, и вектору наблюдения всей 
системы построить математическую модель неизвестной подсистемы. 

Рассмотрим формальную постановку. 

Задан вектор наблюдения 


._ 190 0 
Ч (1) 
где х= р — вектор состояния подсистемы с известной математической моделью, хе"; 


в=[2,...2, | — вектор состояния подсистемы с неизвестной математической моделью, 2е К. 


Математическая модель известной подсистемы — описывается — обыкновенными 
дифференциальными уравнениями: 


х=Р(х, 7, и), (2) 


где и =|и....и ° — векто управления, иеИ< №"; (С — ограниченное множество. 
1 


т 


Для известного управления, действующего в течение заданного интервала времени 
@(+), 1=[0,1, (3) 
определен вектор состояния системы 
г > ах [6 
$(#) =[<(0):#(0]. (4) 
Необходимо найти математическую модель неизвестной подсистемы 
#=5(х, 7, и, Ч), (5) 
[6 
где а= [ч,...4» | — вектор постоянных параметров, че А?; з(х, 2, и, 4) — неизвестная вектор- 
функция. 


о к 0 
Выбор функциональной зависимости #(х, 2, и, 4) и значений параметров = [4..4 | 


должен быть таким, чтобы решение (1) системы совместных дифференциальных уравнений (2) 
и (5) при управлении (3) с начальными условиями 


У(0)=5(0)=[(0):2(0) |", (6) 


минимизировало функционалы 






й „Аа 
(>,-7,) , (7) 


[ 
— 


> = тах | 





У, (1) (1) 

Можно считать, что функциональная зависимость выбирается из счетного множества 
возможных решений. Это означает, что каждой функциональной зависимости должен 
соответствовать целочисленный номер элемента из некоторого множества, который следует 
выбирать, используя методы целочисленной оптимизации. После выбора какой-либо 





Е=ЬИ-Е, 1Е[0,1, ||. (8) 
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функциональной зависимости з(х, 2, и, 4) необходимо определить оптимальные по критериям 


(7) и (8) значения параметров а = т с помощью методов непрерывной оптимизации. Так 


как рассматриваются два критерия оптимизации, следовательно, задача идентификации (1)—(8) 
относится к классу задач многокритериальной смешанной оптимизации. 

В классических методах идентификации [11, 12] функциональную зависимость 
5 (х, 7, и, а) разработчик выбирает на основе анализа практических свойств исследуемой 


системы управления. Сама идентификация предполагает только выбор значений оптимальных 
0 м 
параметров = [4..4 | по методам нелинейного программирования. В последнее время для 


решения подобных задач идентификации используются нейронные сети [13—15], которые 
устанавливаются вместо неизвестной подсистемы (5) и настраиваются на основе обучающих пар 
(3) и (4) по критериям (7) и (8). Решение задачи идентификации на основе нейронной сети 
заведомо предполагает отказ разработчика от попытки найти реальную физическую 
закономерность, которую описывает неизвестная подсистема (5). 

С помощью специальной структуры данных — сетевого оператора — строится множество 


функциональных зависимостей (х, #, и, 4): 


в=|'(х, 2, и, 9): 1=ЬМ}. (9) 
Множество (9) конечно. Для практических задач с помощью сетевого оператора несложно 
построить множество (9) мощностью М> 10. Далее, используя генетический алгоритм, 


. ь : ‚ 0 
выбираем элемент р’ (х, 7, и, 4) и находим оптимальные значения параметров 4’ =[ч1...4, | 5 


Так как в задаче рассматриваем два критерия выбора функциональной зависимости и 
параметров, то решением задачи считаем построение множества Парето в пространстве 
функционалов (7) и (8). Рассмотрим подробнее метод решения. 

Сетевой оператор. Для построения множества функциональных зависимостей (9) введем в 
рассмотрение несколько конечных упорядоченных множеств, из элементов которых состоит 
функциональная зависимость [6, 7]. 


Множество переменных — упорядоченное множество компонент вектора состояний 
объекта кана УЕЁ', 1=1/М№. Множество параметров — упорядоченное множество 
компонент вектора параметров №-= (а, р в а; Е К, 1=Ьр. Множество унарных операций — 


упорядоченное множество функций или однозначных отображений, заданных на числовом 
множестве О, =(р,(2), р (2),....Ри(2)), где р;(2): В! К, У2еВ', ЗуЕВ' > у=р, (2), 


1=ЬЙ’. 
Множество бинарных операций — упорядоченное множество функций двух аргументов 
или однозначных отображений декартового произведения пары одинаковых числовых множеств в 





одно такое же числовое — множество О где 





Х:(2,2"): ЮхВ' = В? >В, \2',2"Е К, ЗуЕА' > у=х, (2. 2"), 1=0,7-1. 

Сетевой оператор — это ориентированный граф, обладающий следующими свойствами: 

а) в графе отсутствуют циклы; 

6) к любому узлу, который не является источником, есть хотя бы один путь от узла- 
источника; 

в) от любого узла, который не является стоком, есть хотя бы один путь до узла-стока; 
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Г) каждому узлу-источнику соответствует элемент из множества переменных У или из 
множества параметров С; 


д) каждому узлу соответствует бинарная операция из множества О, бинарных операций; 
е) каждой дуге графа соответствует унарная операция из множества О; унарных операций. 


Вычисления по сетевому оператору включают три этапа. 

На первом этапе находим узел, который имеет выходящие дуги, но не имеет ни одной 
входящей дуги. Для этой цели берем любой узел и перемещаемся от него к следующему по одной 
из входящих дуг. Продолжаем перемещение до тех пор, пока не найдем узел, у которого нет 
входящих дуг. Первоначально таким узлом является узел-источник. На втором этапе выполняем 
унарную и бинарную операции. Унарная операция соответствует дуге, выходящей из найденного 
узла. В качестве аргумента унарной операции используем значение, которое хранится в 
найденном узле. Бинарная операция соответствует узлу, в который дуга входит. В качестве 
первого аргумента бинарной операции используем либо единицу, либо результат последнего 
вычисления, который хранится в данном узле. В качестве второго аргумента используем 
результат вычисления унарной операции. На третьем этапе исключаем узел и дугу из графа. 
Найденный узел исключаем в том случае, если этот узел не является узлом-стоком и не имеет ни 
одной выходящей дуги. Дугу исключаем во всех случаях, если операция, которая соответствует 
данной дуге, выполнена. 

Этапы вычислений выполняем до тех пор, пока в сетевом операторе не останутся только 
узлы-стоки, в которых хранятся результаты вычислений. 

При вычислении следует соблюдать правило, при котором выполнение унарной операции 
согласно выходящей из узла дуге, осуществляется только в том случае, если в узел, из которого 
данная дуга выходит, не входит ни одна дуга. 

Рассмотрим пример. Пусть задан сетевой оператор (рис. 1). В узлах-источниках расположены 


значения из множества переменных У =(г) и множества параметров № = (<, В, ®, Ф).. В остальных 
узлах расположены номера унарных операций О, = (х ее, 27") ‚ а на узлах 
расположены цифры, соответствующие номерам унарных операций О, (р (2)=2,р.(2)=9щ2, 


рз (2) ==, ра (=)= г") . Тогда, используя правила вычисления результата по сетевому оператору, 


получаем математическое выражение: 


у=ае "зт(0+). 


Рис. 1. Пример сетевого оператора 
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Матрица сетевого оператора — это целочисленная матрица, на диагонали которой 
расположены номера бинарных операций, а остальные элементы либо нули, либо номера 
унарных операций, причем при замене диагональных элементов на нули, а ненулевых 
недиагональных элементов на единицы получаем матрицу смежности графа сети, 
удовлетворяющую условиям «а—в» определения сетевого оператора. 


Для сетевого оператора (см. рис. 1) матрица имеет вид: 


ооооотогоо0 
0000000000 1 
оооо0оо0о00гоо00 
00000100000 
00000010000 
9=0 0000100000. 
00000000020 
00000001300 
00000000040 
00000000011 
00000000001 





Для вычислений результата математического выражения по матрице сетевого оператора 
необходимо определить дополнительные целочисленные векторы: 


— вектор номеров узлов для входных переменных Ъ =[В...Б, |’, где , — номер узла- 





источника в сетевом операторе, с которым связана переменная х,, 1=1, и; 

— вектор номеров узлов для параметров $ = [5..5], где 5, — номер узла-источника в 
сетевом операторе, с которым связан параметр а,, #= Ее 

— вектор номеров узлов для выходных переменных 94 = [4...а„|, где 4, — номер узла 
сетевого оператора, который соответствует выходной переменной у,, 1 =, т. Условие т>1 
соответствует тому, что сетевой оператор описывает не одну формулу. 


Введем вектор узлов а |, размерность Г которого равна количеству узлов 


сетевого оператора. 
Теорема. Пусть задана формула, состоящая из бинарных и унарных операций, и для формул 
построен сетевой оператор, который описывается следующими элементами: матрицей сетевого 





0 
' 


оператора 4 =[ч, |, =, вектором номеров узлов для входных переменных Ъ =[В....Б, | 
вектором номеров узлов для параметров $ =[5...з,|’ и вектором номеров узлов для выходных 
переменных 4 = [4,...4„ |”. Тогда этих элементов достаточно для вычисления формулы. 
Доказательство. Зададим начальные значения вектора узлов 
у, 888 1=Ь,, К=Ьм, 
ап6е 1=5,, /=1., 
„> 8866 1% О ыы а д (е,„ 2) =2, 


где 1=1, [; е, — единичный элемент для бинарной операции ху, . 
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Просматриваем элементы матрицы Ч, находящиеся над главной диагональю. Если 


уу 70, то ‚. Ху, И Ру, (4) п 





Покажем, что 27), 1=4,, К=|)т содержит значения вычислений по формулам, 


которые описывает сетевой оператор. Для этой цели необходимо показать, что выполнены все 
операции графических записей и не нарушен порядок раскрытия скобок. 
Для ненулевого элемента \у, матрицы Ч сетевого оператора выполняются унарная 


операция, соответствующая дуге (1,7), и бинарная операция, соответствующая узлу /, 


следовательно, все операции сетевого оператора будут выполнены, кроме узла 1. Согласно 
топологической сортировке этот узел может быть только источником и с ним не связана бинарная 
операция. 

Пусть в формуле аргументами какой-либо бинарной операции являются унарные операции 


ж&(Рь(7'), р, (2")). В соответствии с топологической сортировкой номер узла /, с которым 
связана бинарная операция А =\у „, должен быть больше номеров узлов { и К, от которых идут 
дуги к узлу / соответствующих унарных операций р =\у;, Г=\р, Поэтому у>Ги ] >К. Таким 
образом, унарные операции будут выполнены раньше, чем бинарная. 

Пусть в формуле аргументом унарной операции является бинарная операция 
р» 6 (2, 2”) К=у,, Р=\,, 1<.7. Следовательно, бинарная операция будет выполнена 


раньше унарной. Согласно алгоритму бинарная операция Хх, (2,2") будет выполняться для всех 


ненулевых элементов \у,, = 0 матрицы сетевого оператора, находящихся в столбце 1, К<р. 
Переход к следующей строке возможен только, если выполнены все операции в 
вышерасположенных строках. Операция р, (х, (2, # будет выполнена при условии, если 
7 


вычислены все унарные операции р. их. (=, =") К <, которые соответствуют всем входящим 


в узел г дугам. Следовательно, при выполнении алгоритма 1 соблюдается правило 
последовательности вычисления унарных операций, порядок раскрытия скобок не нарушается. 
Принцип базисного решения. Для построения алгоритма идентификации системы управления 
используем принцип базисного решения [6, 7]. Принцип заключается в том, что при решении 
оптимизационных задач первоначально строим базисное решение, которое является одним из 
допустимых решений, затем определяем малые вариации базисного решения и создаем алгоритм 
поиска, который ищет оптимальное решение на множестве малых вариаций. 

Определим следующие малые вариации сетевого оператора: 0 — замена унарной 
операции на дуге графа; 1 — замена бинарной операции в узле графа; 2 — добавление в граф 
дуги с унарной операцией; 3 — удаление унарной операции с дугой графа. 

Для описания любой вариации достаточно использовать целочисленный вектор из 


четырех компонент \ =| м м, и; и" ее Первая компонента указывает на номер вариации, вторая 
ВЕ 


— на номер узла, из которого выходит дуга, третья — на номер узла, в который входит дуга, 
четвертая — на номер унарной или бинарной операции. Все вариации сетевого оператора можно 


также выполнять на матрице сетевого оператора 4 =[ч, |, р] =1 Г. Тогда вторая и третья 


компоненты вектора вариаций указывают соответственно на номера строки и столбца матрицы 
сетевого оператора. 
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Пусть задана матрица сетевого оператора 4” = [№4 | ‚ 5/=ЬЁ и вектор вариаций 
[6] 
уу! = Ш ид Ул м | . В результате действия вектора вариаций \' на матрицу сетевого оператора 


4'° получаем новую матрицу сетевого оператора Ч" = [№ | КБЕЬЕ, Ч = оо, 
Каждая вариация изменяет один элемент матрицы 4’ сетевого оператора. Если заданы 
две матрицы сетевых операторов = [№4 |, ЕЕ и + =[4, |, р] =1, Г одинаковой 


размерности Г, то количество несовпадающих между ними элементов определяет количество 
вариаций, которые необходимо выполнить, чтобы из одной матрицы получить другую. 

Для двух несовпадающих матриц сетевых операторов всегда можно найти вариации, 
применение которых к одной матрице позволяет получить другую, причем количество вариаций 
равно расстоянию между матрицами Ф = у“ о...о\! о $9. 

Генетический алгоритм. Для поиска решения задачи идентификации (1)—(8) используем 
генетический алгоритм, построенный на основе принципа базисного решения. 

Первоначально задаем базисное решение задачи, которое может представлять собой 
линейную систему по наблюдаемым координатам пространства состояний или любую другую 
систему, полученную из опыта разработчика системы управления. Базисное решение описываем 


матрицей сетевого оператора 4° = [№ | КБЛ-ЬЕ. 


В качестве возможного решения используем упорядоченное множество векторов вариаций 
и битовую строку, описывающую в коде Грея вектор параметров искомого математического 
ь И И у у ь 0 о 
выражения И” = (м м й 5' аа ‚ =. Н , где Н — количество возможных 
решений или размерность популяции; р — количество параметров в искомом математическом 
выражении; с — число бит целой части; 4 — число бит дробной части. 


Для перевода значений параметров из кода Грея $’ используем следующие соотношения: 


‚ |5, аве (1-1) тоа (с+а)=0, 
где В =”, (9 ) 
7 59, ма. 
Здесь /=1р(с+а). 
Для каждой хромосомы (и, 5 ре 1<1< Н определяем значения функционалов. Для этой 
цели по структурной части Им". ) хромосомы и базисной матрице ЧФ” строим 


матрицу сетевого оператора 4" = уу" о...о\ оФ0. 

Рассматриваемая задача идентификации использует два функционала (7) и (8), поэтому 
решением задачи является множество Парето, для построения которого в генетическом 
алгоритме применяем ранговый критерий. Для каждого возможного решения в популяции 
вычисляем значение ранга доминирования 


Л' = УЛ, (10) 


0 
Здесь У^ и и 
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Решения, имеющие № =0, принадлежат множеству Парето 


Р. = (и. $). 1<А<Н: №=0}. Для репродукции хромосом в популяции случайно отбираем 
пару (и. ") и (и, 5) и определяем вероятность скрещивания по условию 


1+4” 1+4” 


т 11 
1+ А” 1+ 4% (о 


р, = шах 


Где уе [0, 1. 


При скрещивании получаем четыре новых возможных решения (и в =) : 
(и, в), а “”"), (и И В решениях (и. и") ; [ит 8"? сохраняется 
структура родителей, но изменяется параметрическая часть. В решениях (и в ), 


(и в >. изменяются и структурная и параметрическая части. 


Для новых хромосом выполняем операцию мутации с заданной вероятностью ра. 


Алгоритм завершает свою работу после выполнения заданного числа поколений. Результатом 
работы является условное множество Парето, построенное на конечной популяции возможных 
решений. Алгоритм также может завершить свою работу, если условное множество Парето не 
изменяется на протяжении заданного количества поколений. 
Пример. Рассмотрим нелинейную систему управления четвертого порядка. Известная и 
неизвестная подсистемы в математической модели системы имеют порядок пи=А=2. Модель 
известной подсистемы представлена в виде 
ж =х,, 
Хх = -хАчи. 
Наблюдается вектор состояния У=[х Хх. 2) 2. на интервале времени {, =10 с. Вектор 


управления имеет две компоненты и = [м Ш] Идентификация проводится в 21 точке. 


Входной вектор управления имеет следующий вид: 


. 1, 1, 
и ЕЗШ| |, и, =60$| —# |. 
2 2 


При моделировании были выбраны следующие параметры генетического алгоритма: 
размерность популяции — 1024, число поколений — 64, число скрещиваемых пар в одном поколении 
— 256, длина хромосомы — 8, число бит под целую часть — 2, число бит под дробную часть — 14, 
размерность матрицы сетевого оператора 20, число элитарных хромосом — 16, вероятность мутации 
— 0,8, параметр скрещивания — 0,4. 

Заданы следующие конструктивные множества: переменных У = (5 Х, 21, 22, И, и>) р 


ь 1 1-е" 
параметров С=(4,4›), унарных операций д=[н-ь: В 





2? 


— , $12, С0$ 2, агсо2, | и 
1+е 
бинарных операций О, =(2'+ 2", 2'2"). 
В результате получена следующая модель неизвестной подсистемы управления: 
2 =4>2>, 
2, =4121х, 6052) + И> +Х, — 4121 — 422», 


где 4, =1,2091064453125, 4, =0,6346435546875. 
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Для данной модели Л =2,983, Л, =1,0004. 
На рис. 2 приведены результаты моделирования системы управления с полученной 
моделью неизвестной подсистемы. Точками обозначены наблюдаемые значения переменных. 





а) 6) 
|-1 
во 8 
|] о? 
т [|] 
|] Г 
ы о [ ы Г] 
[1 Г 
г [-] 
1 2 3 4 5 6 7 8 9 1 2 3 4 э 6 7 8 
в) 2) 


Рис. 2. Наблюдаемые и полученные значения 


Заключение. В результате использования метода генетического программирования с сетевым 
оператором удалось идентифицировать полностью правую часть подсистемы системы 
управления. Полученная математическая модель представляет собой нелинейные уравнения. 
Решение всей системы удовлетворяет наблюдаемым значениям компонент. 


Работа выполнена по теме гранта РФФИ № 10-08-00618-а. 
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